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Based on a data set obtained in a dental longitudinal study, con- 
ducted in Flanders (Belgium), the joint time to caries distribution 
of permanent first molars was modeled as a function of covariates. 
This involves an analysis of multivariate continuous doubly-interval- 
censored data since: (i) the emergence time of a tooth and the time it 
experiences caries were recorded yearly, and (ii) events on teeth of the 
same child are dependent. To model the joint distribution of the emer- 
gence times and the times to caries, we propose a dependent Bayesian 
semiparametric model. A major feature of the proposed approach is 
that survival curves can be estimated without imposing assumptions 
such as proportional hazards, additive hazards, proportional odds or 
accelerated failure time. 

1. Introduction. The past three decades have witnessed a dramatic de- 
cline in the prevalence of dental caries in children in countries of the Western 
World [De Vos and Vanobbergen (2006)]. However, the disease has now be- 
come concentrated in a small group of children, with the majority unaffected; 
about 10-15% of the children now experience 50% of all caries lesions and 
25-30% suffer 75% of lesions [Marthaler, O'Muhane and Vrbic (1996); Pe- 
tersson and Bratthall (1996)]. The most likely explanation for the difference 
in oral health seems to be socio-economic environmental factors and it oc- 
curs early in childhood [Willems et al. (2005)]. Therefore, to improve dental 
health, early identification of groups at a particular risk of developing caries 
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becomes essential. In this paper we present a Bayesian analysis of a longitu- 
dinal data set, gathered in the Signal- Tandmobiel® study, to investigate the 
relationship between some potential exposure variables and the emergence 
and development of caries in permanent teeth. 

The Signal- Tandmobiel® study is a 6-year longitudinal oral health study 
involving children from Flanders (Belgium) and conducted between 1996 
and 2001. Dental data were collected on gingival condition, dental trauma, 
tooth decay, presence of restorations, missing teeth, stage of tooth eruption, 
orthodontic treatment need, etc. Additionally, information on oral hygiene 
and dietary behavior was collected from a questionnaire completed by the 
parents. The children were examined annually during their primary school 
time by one of sixteen trained and half yearly calibrated dental examiners. 
More details on the Signal- Tandmobiel® study can be found in Section 4.1 
and in Vanobbergen et al. (2000). A primary objective of the investigation 
is to assess the association of some covariates with the emergence and de- 
velopment of caries in permanent teeth. In particular, we are interested in 
studying the effect of the age at start brushing (in years) and of deciduous 
second molars health status [sound/affected; teeth 55, 65, 75, 85, respec- 
tively, see Figure la] on caries susceptibility of the adjacent permanent first 
molars [teeth number 16, 26, 36, 46, see Figure lb]. Additionally, we consid- 
ered the impact of gender (girl/boy), presence of sealants in pits and fissures 
of the permanent first molar (none/present), occlusal plaque accumulation 
on the permanent first molar (none /in pits and fissures / on total surface) and 
reported oral brushing habits (not daily/daily). Note that pits and fissures 




Fig. 1. European notation for the position of (a) deciduous (primary); and (b) per- 
manent teeth. Maxilla = upper jaw, mandible = lower jaw. In (a) the fifth and the eight 
quadrants are at the right-hand side of the subject, and the sixth and the seventh quadrants 
are to the left. In (b) the first and the fourth quadrants are at the right-hand side of the 
subject, and the second and the third quadrants are to the right. 
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Fig. 2. yln example of doubly interval censoring. A scheme of a doubly-mterval-cen- 
sored observation obtained by performing examinations to check the event status at times 
Sill, ■ ■ ■ ,Sii(i. The onset time is left-censored at time ufi = sm, that is, interval- censored 
in the interval {uf'i,uYi] = {0,Siii], the failure time is interval-censored in the interval 



sealing is a preventive action which is expected to protect the tooth against 
caries development. The information on occlusal plaque accumulation, pres- 
ence of sealants in pits and fissures and reported oral brushing habits was 
obtained at the examination where the presence of the permanent first molar 
was first recorded. 

The response of interest is the time to caries development on the perma- 
nent dentition which corresponds to the time from tooth emergence to onset 
of caries. Due to the setup of the study (annual visits of dentists), the onset 
time and the failure time could only be recorded at regular intervals and 
observations on both events were, therefore, interval-censored. A graphical 
illustration of a possible evolution of a tooth is shown in Figure 2. This type 
of data structure, often referred to as doubly-interval-censored failure time 
data, is common in medical research, especially in the context of the analysis 
of acquired immunodeficiency syndrome (AIDS) incubation time, the time 
between the human immunodeficiency virus infection and the diagnosis of 
AIDS. 

Several approaches have been proposed over the past few years for the 
analysis of doubly-interval-censored data. De Gruttola and Lagakos (1989) 
suggested a nonparametric maximum likelihood (NPML) estimator of uni- 
variate survival functions. Alternative methods were subsequently given by 
Bacchetti and Jewell (1991), Gomez and Lagakos (1994), Sun (1995) and 
Gomez and Calle (1999). Kim, De Gruttola and Lagakos (1993) generalized 
the one-sample estimation procedure of De Gruttola and Lagakos (1989) 
to a Cox proportional hazards (PH) model. Their method, however, needs 
to discretize the data. Cox regression with the onset time interval-censored 
and the event time right-censored has been considered by Goggins, Finkel- 
stein and Zaslavsky (1999), Sun, Liao and Pagano (1995) and Pan (2001). 
To simplify the analysis, all of these methods make a rather unrealistic in- 
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dependence assumption between the onset and time-to-event variables [see, 
e.g., Sun, Lim and Zhao (2004)]. 

For the analysis of multivariate doubly-interval-censored survival data, 
frailty models were discussed in Komarek et al. (2005) and Komarek and 
Lesaffre (2008) considering versions of the Cox PH and accelerated fail- 
ure time (AFT) models, respectively. In the latter case, each distributional 
part is specified in a flexible way as a penalized Gaussian mixture with an 
overspecified number of mixture components and under the assumption of 
independence between the onset and time-to-event variables. These models 
provide useful summary information in the absence of estimates of a base- 
line survival distribution and may be formulated in a parametric or semi- 
parametric fashion. However, under these models the regression coefficients 
describe changes in individual responses due to changes in covariates, they 
induce a particular association structure for the clustered variables, and rely 
heavily on the (conditional or subject-specific) assumptions of PH or AFT 
in the relationship between the covariates and the survival times. While the 
PH model assumes the covariates act multiplicatively on a baseline hazard 
function, the AFT model assumes that covariates act multiplicatively on 
arguments of the baseline survival function. Although other type of models, 
such as additive hazards (AH) or proportional odds (PO), could be consid- 
ered in a frailty model context, all these assumptions may be considered 
too strong in many practical applications. For instance, under these models 
survival curves from different covariate groups cannot cross which can be 
unrealistic in some applications [see De lorio et al. (2009)]. This issue is 
particularly relevant for doubly-interval-censored data where the degree of 
available information to perform diagnostic techniques is rather reduced due 
to the censoring mechanism. 

In this paper we discuss a Bayesian semiparametric approach for the anal- 
ysis of multivariate doubly-interval-censored data where the dependence 
across sub-populations, defined by different combinations of the available 
covariates, is introduced without assuming independence between the onset 
and time-to-event variables, without requiring data discretization, and any 
of the commonly used assumptions for the inclusion of covariates in survival 
models. We extend recent developments on dependent nonparametric pri- 
ors, initially proposed by MacEachern (1999, 2000), to provide a framework 
for modeling multivariate doubly-interval-censored data where the resulting 
survival curves have a marginal (or population level) interpretation and are 
not subject-specific. It must be pointed out that the dental data has been 
analyzed before. However, the previous approaches were deficient in that ei- 
ther the doubly-interval-censored nature was not taken into account [Leroy 
et al. (2005)] or restrictive in the sense that the focus was on conditional 
interpretation of the effects of the covariates via frailty models and relying 
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on the AFT or PH assumption [Komarek et al. (2005); Komarek and Lesaf- 
fre (2008)] . Overcoming these problems largely motivates the developments 
presented in this paper. 

The rest of the paper is organized as follows. In Section 2 we introduce 
the proposed model, which is based on the two parameter Poisson-Dirichlet 
process, and discuss its main properties. Section 3 presents the analysis of 
simulated data which illustrate the main advantage of the proposed model. 
Section 4 describes the analysis of the Signal- Tandmobiel® study. A final 
discussion section concludes the article. 

2. The model. 

2.1. Survival regression framework. Let T^j and T.^, i = 1, . . . ,m, j = 
1, . . . , n, be continuous random variables defined on [0, oo) denoting the true 
chronological onset and event times for the jth measurement of the ith 
experimental unit, respectively, and let T--^ = T-^ — T-j be the true time-to- 
event. For example, in our case T^^ is the true time to caries for the jth tooth 
of the ith child, with denoting the true emergence time and T-^ the age 
of caries development. Assume that for each of the m experimental units we 
record the p-dimensional and g-dimensional covariate vectors € X'^ C MP 
and x^- G JY^ C associated to the onset time and to the time-to- 
event , respectively. Let Tf = (I^^, . . . , Tg)' , Tf = (T;f , . . . , TfJ , TJ = 
{Tl...,Tly, T, = (Tf,Tr)', XP = diag(x?;,...,xP„'), Xf = diag(xf;, 
...,x^') and Xi = diag(Xp,X^), i = l,...,m. 

In order to model the joint distribution of the true chronological onset 
times and true time-to-events Tj as a function of covariates, Xj, we consider 

a mixture model. Specifically, we assume Tj|Xj /x^, i = 1, . . . ,m, with 

(2.1) /x.(-|5:,GxJ = j A;2„(V,S)dGx,(/x), 

where k2n{'\f^,'^) denotes a 2n-variate density on with location and 
unstructured scale matrix S taking into account the association among vari- 
ables of the same experimental unit, respectively, and where the mixing 
distributions Gxi , • • • , ^ {Gx : X G X} are dependent probability mea- 
sures. The set of dependent probability measures {Gx : X G X} is defined in 
the complete space of the predictors X and the degree of dependence among 
the elements is governed by the value of the covariates X. If Gx were indexed 
by a finite-dimensional vector of hyper-parameters, for example, normal mo- 
ments, then the model would reduce to a traditional parametric hierarchical 
model. In contrast, in a nonparametric Bayesian approach, every element in 
the set {Gx : X G ^} is a random probability measure and an appropriate 
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prior probability model F for the complete set of unknown distributions in- 
dexed by the set of covariates {Gx : X € X} is specified. In other words, F 
is a distribution over related probability distributions 

(2.2) {Gx:Xe 

Here we focus on the class of discrete random probability measures that can 
be represented as 

oo 

(2.3) G^{B) = Y,^i5eo^)XB), 

1=1 

where i? is a measurable set, wi,W2, • • • are random weights satisfying < 
LOi <1 and -P(X]i^i uji = 1) = 1, and where ^^(x), (') denotes a Dirac measure 
at the random locations 0(X)i,0(X)2, . . . , which are assumed to be inde- 
pendent of the {uji}i^i collection. We discuss specific choices for the random 
probability measure F in (2.2) in the next sections. To better explain our 
proposal, we start with a review of the construction of priors over related 
distributions. 

2.2. Priors over related distributions. The problem of defining priors 
over related random probability distributions has received increasing at- 
tention over the past few years. MacEachern (1999, 2000) proposes the de- 
pendent Dirichlet Process (DDP) as an approach to define a prior model 
for an uncountable set of random measures indexed by a single continuous 
covariate, say, x {Gx :x (z X C M}. The key idea behind the DDP is to cre- 
ate an uncountable set of Dirichlet Processes (DP) [Ferguson (1973)] and to 
introduce dependence by modifying the Sethuraman's (1994) stick-breaking 
representation of each element in the set. If G follows a DP prior with pre- 
cision parameter M and base measure Go, denoted by G '--^ DP{MGo), then 
the stick-breaking representation of G is 

oo 

(2.4) G{B) = Y,^i6e,{B), 

1=1 

where OiIGq '-^^ Go and ui = Vil\j<i{'^ " ^j)^ ^ith Vi\M Beta(l,M). 
MacEachern (1999, 2000) generalizes (2.4) by assuming the point masses 
9{x)i, Z = 1, . . . , to be dependent across different levels of x, but independent 
across I. This approach has been successfully applied to AN OVA [De lorio 
et al. (2004)], survival [De lorio et al. (2009)], spatial modeling [Gelfand, Kot- 
tas and MacEachern (2005)], functional data [Dunson and Herring (2006)], 
time series [Caron et al. (2008)] and discriminant analysis [De la Cruz, Quin- 
tana and Miiller (2007)]. Motivated by regression problems with continuous 
predictors, Grifhn and Steel (2006) and Duan, Guindani and Gelfand (2007) 
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developed models where the dependence is introduced by making the weights 
dependent on covariates. 

Alternatives to these approaches include incorporating dependency by 
means of weighted mixtures of independent random measures [Miiller, Quin- 
tana and Rosner (2004); Dunson and Park (2008)]. This approach was 
originally proposed by Miiller, Quintana and Rosner (2004), motivated for 
the problem of borrowing strength across related submodels. For regression 
problems with continuous predictors, Dunson and Park (2008) proposed a 
countable mixture where the weights depend on the covariates through the 
introduction of a bounded kernel function in the stick-breaking construction 
of the weights. The latter approach requires the choice of a metric for the 
covariate values and, therefore, is not naturally extended to include factors 
and continuous predictors jointly in the model. 

We build our proposal on the construction introduced in De lorio et al. 
(2004) and De lorio et al. (2009) because it is a natural approach to intro- 
duce dependence on both factors and continuous covariates which are com- 
monly of interest in survival models. We consider the class of discrete Linear 
Dependent (LD) models defined as follows. For any given value of the covari- 
ates X G A', in the notation of our motivating problem, the 2n-dimensional 
atoms in the mixing distribution Gx(-) = Z^i^i '^i'^0(X)i (") follow linear (in 
the parameters) models 0(X.)i = Ji.l3i, where the /3;'s represent n{p + q)- 
dimensional vectors of regression coefficients. Therefore, in the dependent 
mixture model given by expression (2.1), P(/i = d(X.)i = X/3;) = ui and the 
dependence is introduced in the point mass locations 0(X)/ through a linear 
model, where the regression coefficients (3i are i.i.d. random vectors from a 

distribution Go, Pi Gq. For simplicity of explanation, consider the case 
of n = 1 and an ANCOVA type of design matrix 



where V is an indicator variable and Z is continuous. For example, V could 
be the gender indicator and Z the age at start brushing. In the LD model 
the dependence across the random distributions is achieved by imposing a 
linear model on the point masses 



As in a standard linear model, (3ii and /S^i can be interpreted as intercepts 
for the point masses associated to the onset time and to the time-to-event, 
respectively, while I32i and /34i are the main effects of gender for the onset and 
time-to-event, respectively, and f3^i can be interpreted as a slope coefficient 
associated to the age at start brushing for the time-to-event. Note that 



X = 



( 



1 y 

1 V z 



) 
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the linear specification is highly flexible and can include standard nonlinear 
transformations of the continuous predictors, for example, additive models 
based on B-splines [see, e.g., Lang and Brezger (2004)], as well as linear 
forms in the continuous predictors themselves. 

2.3. The proposal. In this paper we extend the DDP framework to a 
construction that is based on the general class of Poisson-Dirichlet (PD) 
processes [see, e.g.. Pitman (1996) and Pitman and Yor (1997)]. The PD 
processes belong to the class of species sampling models [see, e.g.. Pitman 
(1996)] and admit the DP prior as an important special case. The PD process 
can also be defined as in expression (4), where the random weights uji are 
independent for the 0;'s and the 6i are i.i.d. from a distribution Gq. The 
weights still admit a stick-breaking representation uji = ViY\j^i{l — Vj), but 

in this case Vj Beta(l — a,b + ja), where either a = —k < and b = 
for some k > and ? = 2, 3, . . . , or < a < 1 and b > —a. We restrict 
our attention to the parameter space A = {ia,b) €M^:0<a<l,6> —a} 
because this is large enough to include two important special cases. When 
a = and b = M, Ferguson's DP{MGo) follows. When a = 7, < 7 < 1, and 
6 = 0, the PD{'y,0) yields a measure whose random weights are based on a 
stable law with index 7. The DP and stable law are key processes because 
they represent the canonical measures of the PD process [Pitman and Yor 



It is now straightforward to extend the Linear Dependent framework to 
the PD process assuming a linear model for the atoms of the process. In this 
way we can define a model for related probability distributions of the form 



where LDPD(a,b,Go) refers to a Linear Dependent PD prior, with parame- 
ters a, b, and Gq. An appealing property of the LDPD survival model given 
by expressions (2.1) and (2.5) is that it can be understood on the basis of 
an equivalent model reformulation as a mixture of multivariate AFT regres- 
sion models. Given a particular matrix of covariates X € A", the vector of 
kernel locations fi in the mixture model (2.1) takes the value X/3, where 
the mixture is defined with respect to the regression coefficients (3. In other 
words, the model can be alternatively formulated by defining the mixture of 
multivariate regression models. 



(1997)]. 



(2.5) 



{Gx : X G X}\a, b, Gq - LDPD{a, b, Go), 



(2.6) 




for all X € and 



(2.7) 



G\a,b,Go-- PD{a,b, Gq). 
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The discrete nature of the PD reahzations leads to their well-known cluster- 
ing properties. The choice of parameters a and b in the PD process controls 
the clustering structure [Lijoi, Mena and Priinster (2007b)]. Given m ob- 
servations, when = (i.e., a DP) the number of clusters n*{m) is a sum 
of independent indicator variables, which implies n*(m)/logm — ?• 6 almost 
surely and n*{m) is asymptotically normal [Korwar and Hollander (1973)]. 
Under the model with < a < 1 and b > —a the sequence {n*(m)} is an 
inhomogeneous Markov chain such that n*{m)/m^ — )• S almost surely, for a 
random variable S with a continuous density on (0, oo) depending on (a, b) 
[Pitman and Yor (1997)]. The asymptotic behavior of the distribution of the 
number of clusters indicates that a general PD model increases as which 
is much faster than the logarithmic rate of the DP model. In general, values 
of a close to 1 favor the generation of a larger number of clusters. 

Besides the clustering structure implied by the extra a parameter in the 
PD process, its role can be also understood when the distribution of PD 
realizations is applied to a partition of the space of interest. In particular, 
for measurable sets B, Bi and B2, with i?i ni?2 = 0, it follows that [Carlton 
(1999)] 

(2.8) \av{G{B)) = Go{B){l-Go{B)) 
and 

(2.9) Cov{G{Bi),G{B2)) = -Go(i?i)Go(i?2) 

Therefore, the extra a parameter controls the variability and covariance of 
disjoint sets of the PD realizations. When a ^ 1, G is highly concentrated 
around Gq and the covariance between disjoint sets is small. When a = we 
recover the corresponding expressions for the DP. Note that the correlation 
between G{Bi) and G{B2) does not depend on the parameter (a, 5) and, 
therefore, is the same as the one arising from the DP model. 

To date, most practical implementations of PD processes have consid- 
ered the parameters a and b as fixed at user-specified values [see, e.g., Ish- 
waran and James (2001)], fixed at empirical Bayes estimates [see, e.g., Lijoi, 
Mena and Priinster (2007a)], or explored the effect of different combinations 
of fixed values for these parameters on the inferences [see, e.g., Navarrete, 
Quintana and Miiller (2008)]. Lijoi, Mena and Priinster (2008), on the other 
hand, proposed independent discrete uniform priors with support points 
{0.01, 0.02, . . . , 0.99} and {0, 1, ... , 2000} for a and 6, respectively. Here we 
allow a and b to be random, having continuous random probability distri- 
butions supported on the restricted parameter space under consideration. 
Moreover, we allow a to be zero with positive probability in order to test 
whether the data arose from LDDP versus a more general LDPD process 
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using a Bayes factor. This additional flexibility can be incorporated at es- 
sentially no additional computational cost. 

2.4. The hierarchical representation. So far, we have focused on mod- 
eling the joint distribution of the survival times of interest, namely, the 
true chronological onset times and true times-to-event T^^ . However, in 
our setting the observed data are given by the events {T^ G (u^-,n^] :i = 
1, . . .,m,j = 1, . . . and {T^^ G (^ij'^ij] = 1> • • • = 1, . . . ,n}, where 
ufj and vfj, and ufj and vfj, represent the lower and upper limits of the 
intervals where the chronological onset, , and event time, , for obser- 
vation j from experimental unit i were observed, respectively. Under the 
assumption of noninformative censoring, we define a model for the events 
Af = {T^^ e {nf^,u^:j = 1, . . . , n} and Af = {i;f G K^.r;,^] : j = 1, • . . ,n}, 
by introducing latent vectors Tf and Tf. We assume 

(2.10) (T?,Tf)|^x/-~-/ix., 

with /ix.(Tf,Tf|S,G) = /x,(Tp,Tf-TP|5],G) and where /x,(-|S, G) is 
defined as in (2.6). Notice that a choice of the continuous kernel k defines 
the model. A multivariate log-normal distribution is convenient for prac- 
tical reasons. Let Zj = (logT/^, . . . ,logr.^,logr-^, . . . , logT-^)' denote the 
logarithmic transformation of the true chronological onset times and true 
times-to-event such that 

(2.11) /x,(Ti|5],G) = I ^N2n{zi I X,/3,S)n^,7i^ dG{f3), 

where A'2n('|/^, S) refers to a 2n-dimensional normal distribution with mean 
H and covariance matrix H. The mixture model /x^ can be equivalently 
written as a hierarchical model by introducing latent variables /3* such that 



(2.12) Zi|/3*,S'-^''-Af2n(X,/3*,S), 

(2.13) /3t,...,/3:;|G'~-G 
and 

(2.14) G\a,b,Go-- PD{a,b,Go), 



where the baseline distribution Go is assumed to be n(p + g)-dimensional 
normal distribution Go(/3) = Nn(p-\-q){ia, S). 
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2.5. Some properties. An important property of the proposed model 
given by expressions (2.11)-(2.14) is that the complete distribution of sur- 
vival times is allowed to change with values of the predictors (including 
properties such as skewness, multimodality, quantiles, etc.) instead of just 
one or two characteristics, as implied for many commonly used survival 
models. However, we make explicit the dependence of some functionals of 
interest of the distribution of the event times on the covariates in order to 
compare them to the corresponding expression arising from the commonly 
used models. The implied marginal mean, hazard function and cumulative 
distribution (CDF) function for coordinate j in the vector Tj, Tjj, as func- 
tions of the associated vector of the design matrix Xj , Xjj , are given by 

oo 

(2.15) ^(r,,|xy) = J^a;iexp{x^A + 0.5a2}, 

1=1 



(2.16) H|x.W 
and 



E£i '^i/o.al (exp{-x^/3 Jt) 



OO 



(2.17) ^T,,|x,,(i) = 5^wzFo_,2(exp{-x^/3Jt), 

1=1 

respectively, where /q ^2 and Fq ^i refers to the density and CDF of a lognor- 
mal distribution with mean and variance o"^ , and (t| = S jj . These expres- 
sions show the additional flexibility associated to the proposed model. For 
instance, in contrast to a simple AFT survival model based on the lognormal 
distribution, the mean function of our proposal given by expression (2.15) 
is a convex combination of exponential functions. Furthermore, the implied 
CDF given by expression (2.17) is a convex combination of CDF's arising un- 
der the AFT model, FT..^^..{t) = Fq^2 (exp{— x^j/3}t), where covariates act 

multiplicatively on arguments of the baseline survival function. This simple 
fact induces an important property of our proposal, namely, that survival 
curves are allowed to cross for different values of a predictor, which is not 
possible under the AFT assumption. Other commonly used models such as 
PH, AH and PO will also fail to capture this behavior. Under the PH, AH 
and PO models, the dependence of the CDF on predictors is given by 

l-F^^^|.^^(t) = {l-Fo,,.(t)rPK^>, 

1 - F7^,^ix,,(i) = {1 - Fo_,2(t)}exp{-x^/3t} 

and 



■exp{x/3}. 
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respectively. Notice that this constraint associated to the commonly used 
models remains if i^o.o-j is modeled in a nonparametric manner and/or if the 
linear form is replaced for a more general function m(xjj). Although 
some fixes have been proposed in the context of PH models for this un- 
appealing property, for example, the inclusion of interactions with time or 
stratification, our modeling approach has proved to be a more flexible alter- 
native. We refer to De lorio et al. (2009) for a thorough comparison in the 
context of univariate (not doubly censored) survival data. 

2.6. Prior distributions and MCMC implementation. For a and h we 
consider joint prior distributions of the kind p{a, b) = p{a)p{b\a), where p{a) 
is a mixture of point mass at zero and a continuous distribution on the unit 
interval (0, 1) and p{b\a) is a continuous distribution supported on {—a,oo). 
More specifically, we assume 

(2.18) a|A, ao, ai ~ A5o(") + (1 ~ '^) Beta(-|ao, ai) 
and 

(2.19) b\a,fMb,crb ~ N{fMb,crb)I{-a, oo), 

where < A < 1, and Beta(-|ao, ai) refers to a beta distribution with pa- 
rameters ao and ai. This modeling strategy allows us to explicitly compare 
a DP model versus an encompassing PD alternative. Notice that this is an 
important component because the evaluation of any other model comparison 
criteria would require the computation of a highly complex area under the 
multivariate normal distribution which is difficult to be performed in prac- 
tice. Finally, to complete the model specification, we assume independent 
hyper-priors m ~ iV„(p+g)(?7, T), S ~ /W„(p+,)(7,r), and S ~ /W2n(j^, ^^), 
where IW2ni'^,^) denotes a 2n-dimensional inverted- Wishart distribution 
with degrees of freedom and scale matrix fi. 

The hierarchical representation of the model allows straightforward pos- 
terior inference with Markov Chain Monte Carlo (MCMC) simulation. As in 
the context of standard DP models, two different kinds of MCMC strategies 
could be considered for computation in the LDPD model: (I) to marginalize 
out the unknown infinite-dimensional distributions [see, e.g., Ishwaran and 
James (2003) and Navarrete, Quintana and Miiller (2008)] or (II) to employ 
a truncation to the stick-breaking representation of the process [see, e.g., 
Ishwaran and James (2001)]. In the case (I), several alternative algorithms 
could be considered to sample the cluster configurations: (La) via a Gibbs 
scheme through the coordinates [see Navarrete, Quintana and Miiller (2008) 
for a discussion in the PD context] or (I.b) to adapt reversible-jump- like 
algorithms [see, e.g., Dahl (2005)] to the PD context. Functions implement- 
ing these approaches were written in a compiled language and incorporated 
into the R library "DPpackage" [Jara (2007)]. A complete description of the 
full conditionals and algorithms is available in the supplemental article [Jara 
et al. (2010a)]. 
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3. An illustration using simulated data. To validate our approach, we 
conducted the analysis of real-life and simulated data sets. The results of the 
real-life data analysis are reported in the supplemental article [Jara et al. 
(2010b)]. The simulated data sets mimic to a certain extent the Signal- 
Tandmobiel® data. We consider one onset time Tf and one time-to-event 
time for m = 500 subjects. We assume a binary predictor and 250 sub- 
jects in each level (groups A and B). Different distributions were assumed 
for each level of the predictor such that 



log(T«,Tf),. 



O rpT 

250 ! 250 



/a 



and 



log(T, 



O rpT 

251) -^251 



rpT 

' 500^^ 500 



,log(T5^ 

Two scenarios for the distributional parts of the model were considered. In 
scenario 1, a mixture of two bivariate lognormal distributions was assumed 
for group A while a bivariate lognormal distribution was assumed for group 
B. An important characteristic of scenario 1 is the bimodal behavior of the 
distribution of the onset time and time-to-event in group A. In group B, a 
unimodal behavior for the distribution of both variables was assumed. In 
scenario 11, mixtures of bivariate lognormal distributions were assumed for 
both groups. However, the components of the mixtures were specified in 
such a way that, for group A, the onset times follow a bimodal distribu- 
tion and the time-to-events follow a unimodal distribution. In group B, the 
reverse behavior was assumed, namely, the onset times follow a unimodal 
distribution while the time-to-events a bimodal distribution. 

In both scenarios and variables of interest, the survival curves for both 
groups cross. The true distributions in each scenario are given next: 

• Scenario I: Mixture model for group A-Single model for group B. 
fA = 0.5 X N2 



+ 0.5 X iV2 





"1.80' 


,10"^ 


"5.00 


2.50" 


( 


0.75 


2.50 


300 



2.40 
3.00 



,10' 



2.50 
1.25 



1.25 
100 





"2.1" 


,10~2 


"3.24 


8.10" 


( 


2.2 


8.10 


64 



and 

fB = N2 

Scenario II: Mixture model for both groups A and B 
fA = 0.5 X N2 

+ 0.5 X N2 





"1.8" 


,10"^ 


"5.50 


2.50" 


( 


2.2 


2.50 


640 



2.4 
2.2 



,10- 



2.50 
1.25 



1.25 

640 
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and 



2.10 
0.75 



,10' 



3.24 
8.10 



8.10 
30.00 





"2.10' 


,10"^ 


"32.4 


1.25" 


( 


0.75 


1.25 


100 



Jb ^ 0.5 X N2 



+ 0.5 X N2 



The true onset and event times were interval-censored by simulating the 
visit times for each subject in the data set. The first visit was drawn from 
an A^(7,0.2^) distribution. Each of the distances between the consecutive 
visits was drawn from an A^(l,0.05^) distribution. 

The LDPD model was fitted to both simulated data sets using the follow- 
ing values for the hyper-parameters: A = 0.5, ao = «i = 1) f^b = 10, (Jb = 200, 
u = A, f2 = I2, 7 = 5, r = I4, 77 = O4 and T = IOOI4. In each analysis 4.02 
millions of samples of a Markov chain cycle were completed. Because of 
storage limitations and dependence, the full chain was subsampled every 
200 steps after a burn-in period of 20,000 samples, to give a reduced chain 
of length 20,000. 

Figures 3 and 4 display the true and estimated survival curves for the 
onset and time-to-event under scenarios I and II, respectively. The predictive 
survival function closely approximated the true survival functions, which 
were almost entirely enclosed in pointwise 95% highest posterior density 
(HPD) intervals. We note that these results are for one random sample from 
two particular densities, and these conclusions should not be overinterpreted. 
Nonetheless, these examples do show that our proposal is highly flexible 
and is able to capture different behaviors of the onset and time-to-event 
survival functions. The examples also show that when a parametric model 
is appropriated, the proposed model does not overfit the data. 



4. The Signal- Tandmobiel® data. 

4.1. The Signal- Tandmobiel^ study and the research questions. For this 
project 4468 children were examined on a yearly basis during their primary 
school time (between 7 and 12 years of age) by one of sixteen dental ex- 
aminers. Sampling of the children was done according to a cluster-stratified 
approach with 15 strata. A stratum consists of a particular combination of 
one of the five provinces in Flanders with one of the three school systems. 
Schools were selected such that all children had equal probability of being 
selected and for each school all children of the first class were examined. 
Clinical data were collected by the examiners based on visual and tactile 
observations (no X-rays were taken), and data on oral hygiene and dietary 
habits were obtained through structured questionnaires completed by the 
parents. 
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10 30 50 10 30 50 

Time Time 
(c) (d) 

Fig. 3. Simulated data — Scenario 1; Estimated survival functions for the onset and time- 
to-event times for the group A are displayed m panels (a) and (c), respectively. Estimated 
survival functions for the onset and time-to-event times for the group B are displayed in 
panels (b) and (d), respectively. The posterior means (solid lines) are presented along the 
pointwise 95% HPD intervals. The true functions are presented in dashed lines. 

The primary interest of our analysis is to study the relationship between 
age at start brushing (in years) and deciduous second molars health status 
(sound/affected) with caries susceptibility of the adjacent permanent molars. 
Here, "affected molar" refers to a tooth that is decayed, filled or missing due 
to caries. The deciduous second molars refer to teeth 55, 65, 75 and 85 and 
first molars refer to teeth 16 and 26 on the maxilla (upper quadrants), and 
teeth 36 and 46 on the mandible (lower quadrants). The numbering of the 
teeth follows the FDI (Federation Dentaire Internationale) notation which 
indicates the position of the tooth in the mouth (see Figure 1). Position 26, 
for instance, means that the tooth is in quadrant 2 (upper left quadrant) and 
position 6 where numbering starts from the mid-sagittal plane. The level of 
decay was scored in four levels of lesion severity: dA (dentine caries with 
pulpal involvement), 32 (limited dentine caries), d2 (enamel cavity) and dl 
(white or brown-spot initial lesions without cavitation). Here we consider 
level (is of severity, which defines a progressive disease. 

Note that for about five years the deciduous second molars are in the 
mouth together with the permanent first molars. It is thus possible that a 
caries process on the primary and permanent molar occurs simultaneously. 
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10 30 50 10 30 50 

Time Time 
(c) (d) 

Fig. 4. Simulated data — Scenario 2; Estimated survival functions for the onset and time- 
to-event times for the group A are displayed in panels (a) and (c), respectively. Estimated 
survival functions for the onset and time-to-event times for the group B are displayed in 
panels (b) and (d), respectively. The posterior means (solid lines) are presented along the 
pointwise 95% HPD intervals. The true functions are presented in dashed lines. 

In this case it is difficult to know whether caries on the deciduous molar 
caused caries on the permanent molar or vice versa. For this reason, the 
permanent first molar was excluded from the analysis if caries were present 
when emergence was recorded. Moreover, the permanent first molar had to 
be excluded from the analysis if the adjacent deciduous second molar was 
not present in the mouth already at the first examination. For 948 children 
none of the permanent first molars was included in the analysis due to the 
previously mentioned reasons. In total, 3520 children (12,485 permanent first 
molars) were included in the analysis of which 187 contributed one tooth, 
317 two teeth, 400 three teeth and 2616 all four teeth. 

4.2. The analysis and the results. We consider gender (0 = boy, 1 = girl) 
and the status of the adjacent deciduous second molar (sound = 0, affected = 
1) as covariates for the emergence times Tj*?, namely, to define the design 
vectors x^. For the time-to-caries variables, we use a similar set of covariates 
as Leroy et al. (2005), namely, the covariate vectors x^- for the caries part of 
the model include gender, presence of sealants on the permanent first molar 
(0 = absent, 1 = present), occlusal plaque accumulation for the permanent 
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first molar (0 = none, 1 = in pits and fissures or on total surface), reported 
oral brushing habits (0 = not daily, 1 = daily) and status of the adjacent 
deciduous second molar. In contrast to Leroy et al. (2005), we did not use 
the status of the adjacent deciduous first molar as a covariate due to its 
large dependence on the status of the adjacent deciduous second molar and 
included the age at start brushing in a linear fashion. 

For the model, 4.02 millions of samples of a Markov chain cycle were 
completed. Because of storage limitations and dependence, the full chain 
was sub-sampled every 200 steps after a burn-in period of 20,000 samples, 
to give a reduced chain of length 20,000. We consider A = 0.5 reflecting equal 
prior probabilities for the LDDP and LDPD models. The values of the other 
hyper-parameters were taken as ao = ai = 1, = 10, (Jb = 200, u = 10, 
S7 = Is, 7 = 31, r = I28, ^7 = O28 and T = 100 x I28. We also performed 
the analysis with different hyper-parameters values, obtaining very similar 
results. This suggests robustness to the prior specification. 

The posterior probability for a = was 21.63%. Correspondingly, the 
Bayes factor for the hypothesis of a LDPD against the DP version of the 
model was 3.62. This result suggests a "substantial" support of the data 
to the PD version of the model according to the Jeffreys' scale [Jeffreys 
(1961), page 432]. As Bayes factors may be sensitive to the prior specifica- 
tion, we performed a sensitivity analysis using different prior weights on the 
LDDP versus a more general LDPD model. Specifically, we chose A = 0.3 
and A = 0.7. The corresponding Bayes factors for the LDPD against the DP 
version of the model were 2.72 and 2.21, respectively. The results, therefore, 
indicate robustness of the model choice to the prior specification. More im- 
portantly, in all cases the PD version of the model is to be preferred when 
compared to the single precision DP model. 

The emergence and caries processes showed a nonsignificant association, 
evaluated by the Pearson correlation coefficient on the log-scale induced by 
5], for most of the teeth, except for tooth 46 where a small negative associ- 
ation was observed. The posterior mean (95% HPD intervals) for the emer- 
gence and caries processes for tooth 16, 26, 36 and 46 were —0.06 (—0.18; 
0.05), -0.06 (-0.18; 0.07), -0.05 (-0.13; 0.02) and -0.10 (-0.18; -0.02), 
respectively. The association among emergence times and among time-to- 
caries was positive and significant. Table 1 displays the posterior means and 
95% HPD intervals for the Pearson correlation among the teeth. The re- 
sults indicate an exchangeable correlation matrix would suffice to explain 
the emergence process. However, this type of association structure does not 
hold for the caries process. The Pearson correlation was bigger for the log 
time-to-caries for teeth in the same jaw. Similar and lower associations were 
observed when considering diagonally or vertically opponent teeth. Thus, 
the results suggest that the correlation structure induced for frailty models 
is not appropriate for these data. 
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Table 1 

Signal- Tandmobiel® study: Posterior mean (95% HPD interval) for the Pearson 
correlation coefficient between log emergence times (upper diagonal) and log 
time-to-caries (lower diagonal) for different teeth 



Tooth 




Tooth 




16 


26 


36 


46 


16 




0.60 (0.56; 0.64) 


0.60 (0.56; 0.64) 


0.60 (0.56; 0.64) 


26 


0.88 (0.81; 0.94) 




0.59 (0.55; 0.63) 


0.59 (0.57; 0.63) 


36 


0.47 (0.35; 0.57) 


0.43 (0.30; 0.55) 




0.61 (0.57; 0.65) 


46 


0.44 (0.28; 0.61) 


0.39 (0.22; 0.58) 


0.61 (0.54; 0.67) 





In contrast to NPML approaches, an important characteristic of the pro- 
posed model is the abihty to make inferences on any quantile of interest. 
With respect to the median, neither the emergence nor the caries process 
exhibit a significant difference among the four permanent first molars. For 
all combinations of covariates, molars of girls tend to emerge earlier than 
those of boys. However, nonsignificant differences were found. Regarding 
caries experience, the difference between boys and girls was not significant, 
however, the frequency of brushing, presence of sealant, presence of plaque, 
age at start brushing and caries experience of neighboring deciduous second 
molars have a significant effect on the caries process. Table 2 shows the pos- 
terior mean and the 95% HPD interval for the median emergence time and 
time-to-caries for teeth 36 and 46 of boys with the "best," "worst" and two 
intermediate combinations of discrete covariates. The results are shown for 
4 different values of age at start brushing. 

Figures 5 and 6 illustrate the estimated hazard and survival functions 
for the time-to-caries for tooth 16 in boys with the "best," "worst" and 
two intermediate combinations of the discrete covariates by age at start 
brushing. For children who started brushing their teeth after the age of 5, 
a high peak in the hazard function of caries is observed already less than 
1 year after emergence. A smaller peak, shifted to the right and of much 
lower magnitude, was observed for children who brush their teeth before the 
age of 5. Furthermore, for a given combination of the discrete predictors, 
the hazard function for caries crossed for different values of age at start 
brushing, suggesting that a proportional hazards model is not an appropriate 
alternative for modeling the time to caries. For a given age at start brushing, 
the presence of an affected deciduous second molars significantly increases 
the pick in the hazard function of caries in the permanent first molar. When 
the teeth are daily brushed since an early age, plaque-free and sealed the 
hazard for caries starts to increase approximately 2 years after emergence, 
whereas when the teeth are not brushed daily and are exposed to other risk 



Table 2 

Signal- Tandmobiel® study: Posterior mean (95% HPD interval) for the median emergence time and time-to- caries since emergence 
(years) for some covariate combinations and teeth. The results are shown for boys and teeth 36 and 46 with the following combination of 
the covariates: Gl for no plaque, present sealing, daily brushing and sound primary second molar, G2 for no plaque, present sealing, 
daily brushing and affected primary second molar, G4 for present plaque, no sealing, not daily brushing and sound primary second 
molar, and G4 for for present plaque, no sealing, not daily brushing and affected primary second molar 



Age at start 
brushing (years) 



Emergence 



Caries 



Covariate group 



Tooth 36 



Tooth 46 



Tooth 36 



Tooth 46 



1 


Gl 


6.57 


(6.54 


6.60) 


6.56 


(6.53 


6.60) 


12.62 


(11.44; 13.82) 


11.89 


(10.65; 13.17) 




G2 


6.58 


(6.54 


6.61) 


6.57 


(6.54 


6.61) 


9.99 


(8.80 


11.18) 


9.72 


(8.45 


11.04) 




G3 


6.57 


(6.54 


6.60) 


6.56 


(6.53 


6.60) 


7.72 


(6.68 


8.54) 


8.49 


(6.95 


9.79) 




G4 


6.58 


(6.54 


6.61) 


6.57 


(6.54 


6.61) 


5.98 


(4.98 


6.85) 


6.83 


(5.49 


7.94) 


3 


Gl 


6.57 


(6.54 


6.60) 


6.56 


(6.53 


6.60) 


11.08 


(9.82 


12.29) 


10.48 


(9.24 


11.765) 




G2 


6.58 


(6.54 


6.61) 


6.57 


(6.54 


6.61) 


8.63 


(7.65 


9.73) 


8.47 


(7.23 


9.63) 




G3 


6.57 


(6.54 


6.60) 


6.56 


(6.53 


6.60) 


6.66 


(5.85 


7.46) 


7.37 


(6.32 


8.39) 




G4 


6.58 


(6.54 


6.61) 


6.57 


(6.54 


6.61) 


5.16 


(4.38 


5.94) 


5.94 


(5.04 


6.75) 


5 


Gl 


6.57 


(6.54 


6.60) 


6.56 


(6.53 


6.60) 


9.67 


(8.09 


11.28) 


9.25 


(7.39 


11.29) 




G2 


6.58 


(6.54 


6.61) 


6.57 


(6.54 


6.61) 


7.49 


(6.32 


8.72) 


7.47 


(5.86 


9.18) 




G3 


6.57 


(6.54 


6.60) 


6.56 


(6.53 


6.60) 


5.78 


(4.85 


6.74) 


6.47 


(5.33 


7.65) 




G4 


6.58 


(6.54 


6.61) 


6.57 


(6.54 


6.61) 


4.47 


(3.71 


5.31) 


5.22 


(4.22 


6.20) 


7 


Gl 


6.57 


(6.54 


6.60) 


6.56 


(6.53 


6.60) 


8.46 


(6.50 


10.45) 


8.28 


(5.69 


11.21) 




G2 


6.58 


(6.54 


6.61) 


6.57 


(6.54 


6.61) 


6.54 


(5.07 


8.01) 


6.69 


(4.56 


9.11) 




G3 


6.57 


(6.54 


6.60) 


6.56 


(6.53 


6.60) 


5.04 


(3.91 


6.25) 


5.76 


(4.26 


7.53) 




G4 


6.58 


(6.54 


6.61) 


6.57 


(6.54 


6.61) 


3.91 


(3.00 


4.87) 


4.65 


(3.38 


6.14) 
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to 



Years since emergence 

(c) 




Years since emergence 
(d) 



Fig. 5. Signal- Tandmohiel® study: Estimated hazard function for tooth 16 of hoys who 
started brushing their teeth at the age of 1 (solid line), 3 (dashed line), 5 (dotted line) or 7 
(dotted-dashed line). Panels (a) and (b) present the results for no plaque, present sealing, 
daily brushing and sound primary second molar (a) or affected primary second molar (b) . 
Panels (c) and (d) present the results for present plaque, no sealing, not daily brushing 
and sound primary second molar (c) or affected primary second molar (d) . 



factors the hazard starts to increase immediately after emergence. The peak 
in the hazard for caries after emergence can be explained by the fact that 
teeth are most vulnerable for caries soon after emergence when the enamel 
is not yet fully developed. The curves for girls were similar, and are therefore 
omitted. 

Figure 6 also shows the way in which the age at start brushing is related 
to the caries process. The bigger the age at start brushing, the bigger the 
prevalence of caries. However, this increase in the prevalence is only ob- 
served in the first years after emergence. After 5 years since emergence, the 
prevalence of caries experience tends to be the same (and can in fact be 
the same, depending on the exposure to other risk factors) regardless of the 
age at start brushing. This result suggests that PH, AFT, AH or PO mod- 
els are not appropriate for the analysis of caries experience since their are 
constrained in such a way that survival curves are not allowed to cross for 
different values of a predictor. Although the peak in the hazard for caries at 
approximately 1-2 years after emergence was also observed in Leroy et al. 
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Years since emergence 
(a) 



Years since emergence 
(b) 



Years since emergence 
(c) 



Years since emergence 
(d) 



Fig. 6. Signal- Tandmobiel® study: Estimated survival function for tooth 16 of hoys who 
started brushing their teeth at the age of 1 (solid line), 3 (dashed line), 5 (dotted line) or 7 
(dotted-dashed line). Panels (a) and (b) present the results for no plaque, present sealing, 
daily brushing and sound primary second molar (a) or affected primary second molar (b) . 
Panels (c) and (d) present the results for present plaque, no sealing, not daily brushing 
and sound primary second molar (c) or affected primary second molar (d) . 



(2005) and Komarek and Lesaffre (2008), this interesting finding was not 
detected due to the models considered by these authors. 

5. Concluding remarks. We have introduced a probability model for de- 
pendent random distributions in the context of multivariate doubly-interval- 
censored data. The main features of the proposed model are ease of inter- 
pretation, the ability of testing the hypothesis of the independence between 
onset and time-to-event variables, efficient computation and the fact that 
assumptions on survival curves, such as proportional hazards, additive haz- 
ards, proportional odds or accelerated failure time, are not needed. 

The proposal is based on a LDPD model, which contains the LDDP model 
as an important special case, and is specified in such a way that a simple 
hypothesis test for a LDDP versus a more general LDPD alternative can 
be performed with no real additional computational effort and without the 
need of independent fit of the models. 
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Several extensions of this work are possible. We are currently working on 
a version of the model that takes into account potential misclassification of 
the caries process and its effect on the corresponding inferences. Finally, the 
extension of the model allowing for weight dependent covariates is also the 
subject of ongoing research. 
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SUPPLEMENTARY MATERIAL 

Supplement A: MCMC schemes for posterior computation 

(DOI: 10.1214/10-AOAS368SUPPA; .pdf). A complete description of the 
full conditionals for marginal and conditional MCMC algorithms for fitting 
the LDPD survival model for doubly-interval-censored data is given. 

Supplement B: The HIV-AIDS data (DOI: 10.1214/10-AOAS368SUPPB; 
.pdf). The analysis of the data set considered by De Gruttola and Lagakos 
(1989) is presented. This analysis allows for the comparison of the LDPD 
model with the one-sample nonparametric maximum likelihood estimator 
proposed by De Gruttola and Lagakos (1989). The data set considers infor- 
mation from a cohort of hemophiliacs at risk of human immunodeficiency 
virus (HIV) infection from infusions of blood they received periodically to 
treat their hemophilia in two hospitals in France. For this cohort both in- 
fection with HIV and the onset of acquired immunodeficiency syndrome 
(AIDS) or other clinical symptoms could be subject to censoring. There- 
fore, the induction time between infection and clinical AIDS are treated as 
doubly-censored. 
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